Calculate the RMS between two observatories



[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/image/work/IMAGE_ARRAY.xlsx');
master_b(1,:) = [];

time_array_min = (datenum(1995,1,1,0,0,30): (1/(24*60)): datenum(2012,12,31,23,59,30))';
a = length(time_array_min);
S = dir('/Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/*.mat');

for i = 1:length(S)-1,
    
    eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(i ,1)) '.mat' ' x_data y_data z_data']);
    
    x_data_1 = x_data;
    y_data_1 = y_data;
    z_data_1 = z_data;
    
    
    for j = i+1:length(S),
        
        metadatacellarray = master_b(i ,:);
        lat = master_a(i ,1);
        long = master_a(i ,2);
        alt = master_a(i ,3);
        eval(['load /Users/manojnair/projects/obs_mag_data/image/work/station_matlab_files/' cell2mat(master_b(j ,1)) '.mat' ' x_data y_data z_data']);
        
        
    end;
end;


% steps

% read x,y, z data
% calculate F
% read x,y,z secular removed data
% find scalar F through IGRF
% compare the above two

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);

nc_fname = '/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/ALE_1995_2012_Absolute.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);

% fit spline to F and remove the secular variations

f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);


L = isnan(f_data);

data_array = f_data(~L);

time_array = time_array_min(~L);


b1 = min(time_array):365:max(time_array)+10;

if length(b1) > 1,
    
    sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
    
    v=ppval(time_array_min,sp);
    
else,%data length <=1 year
    
    sp=robustfit(time_array,data_array);
    
    v= sp(1) + sp(2) * time_array_min;
    
end;


f_data_secular_removed = f_data - v';


% now read the secular removed data and calculate F using IGRF

nc_fname = '/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/ALE_1995_2012_Secular_Removed.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);


icount = 1;
for decyear = 1995:0.1:2013,
    [mag_igrf(icount,:)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
    icount = icount + 1;
end;

L = isnan(x_data) | isnan(z_data) | isnan(y_data);

% get values at all points through  interpolation
MainB = interp1(datenum(1995:0.1:2013,1,1), mag_igrf, time_array_min, 'linear');

% get a unit vector in the direction of Main field

Unit_Vector = [MainB(~L,:)./repmat(sqrt(MainB(~L,1).^2+MainB(~L,2).^2+MainB(~L,3).^2),[1,3])];

% find dot product

scalar_f = dot([x_data(~L) y_data(~L) z_data(~L)]',Unit_Vector');


% May 16,2013. Both scalar_f and f_secular_removed are the same



%% load all the secular removed data , find F and save the data.
% While the secular removed f can be calculated from absolute X,Y,Z and
% then spline fitting, this is better (as I found with the above
% experiment)

S = dir('/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/*.nc');
S1 = dir('/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/*.nc');

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);


for i = 1:length(S),
    
    
    nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/' S(i).name] ;
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
    
    
    icount = 1;
    for decyear = 1995:0.1:2013,
        [mag_igrf(icount,:), H(icount), DEC(icount), DIP(icount)] = igrf11magm(0,obj.geospatial_lat,obj.geospatial_lon,decyear);
        icount = icount + 1;
    end;
    
    L = isnan(x_data) | isnan(z_data) | isnan(y_data);
    
    % get values at all points through  interpolation
    MainB = interp1(datenum(1995:0.1:2013,1,1), mag_igrf, time_array_min, 'linear');
    declination = interp1(datenum(1995:0.1:2013,1,1), DEC, time_array_min, 'linear');
    inclination = interp1(datenum(1995:0.1:2013,1,1), DIP, time_array_min, 'linear');
    
    
    % get a unit vector in the direction of Main field
    
    Unit_Vector = [MainB(~L,:)./repmat(sqrt(MainB(~L,1).^2+MainB(~L,2).^2+MainB(~L,3).^2),[1,3])];
    
    % find dot product
    
    scalar_f = nan(size(x_data));
    
    scalar_f(~L) = dot([x_data(~L) y_data(~L) z_data(~L)]',Unit_Vector');
    
    % compute declination
    
    
    declination_est = nan(size(x_data));
    
    L = isnan(x_data) | isnan(y_data);
    
    declination_est(~L) = 180/pi * atan((MainB(~L,2) + y_data(~L))./(MainB(~L,1) + x_data(~L)));
    
    % remove annual variation from the declination data
    
    data_array = declination_est(~L);
    
    
    declination_est_secular_removed = declination_est - declination';
    
    
    % compute inclination
    
    
    inclination_est = nan(size(x_data));
    
    L = isnan(x_data) | isnan(y_data) | isnan(z_data);
    
    
    inclination_est(~L) = 180/pi * atan((MainB(~L,3) + z_data(~L) )./ ...
        sqrt( (MainB(~L,2) + y_data(~L)).^2 + (MainB(~L,1) + x_data(~L)).^2 ) );
    
    
    % remove annual spline from the inclination data
    
    
    inclination_est_secular_removed = inclination_est - inclination';
    
    
    % Now if available calculate the F, declination and inclination in the traditional way
    
    nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Absolute_netCDF/' S(i).name(1:14) 'Absolute.nc'];
    
    
    if (~isempty(dir(nc_fname))),
        
        [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
        
        % fit spline to F and remove the secular variations
        
        f_data = sqrt(x_data.^2 + y_data.^2 + z_data.^2);
        
        
        L = isnan(f_data);
        
        data_array = f_data(~L);
        
        
        
        time_array = time_array_min(~L);
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v=ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp=robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        f_data_secular_removed = f_data - v';
        
        % Calculate declination and remove the baseline
        
        
        declination_est = nan(size(x_data));
        
        L = isnan(x_data) | isnan(y_data);
        
        declination_est(~L) = 180/pi * atan( y_data(~L)./ x_data(~L) );
        
        L = isnan(declination_est);
        
        data_array = declination_est(~L);
        
        
        
        time_array = time_array_min(~L);
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v=ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp=robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        declination_secular_removed = declination_est - v';
        
        
        % Calculate inclination and remove the baseline
        
        
        inclination_est = nan(size(x_data));
        
        L = isnan(x_data) | isnan(y_data) | isnan(z_data);
        
        
        inclination_est(~L) = 180/pi * atan( z_data(~L) ./ ...
            sqrt(  y_data(~L) .^2  + x_data(~L).^2  ) );
        
        
        
        L = isnan(inclination_est);
        
        data_array = inclination_est(~L);
        
        
        
        time_array = time_array_min(~L);
        
        
        b1 = min(time_array):365:max(time_array)+10;
        
        if length(b1) > 1,
            
            sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
            
            v=ppval(time_array_min,sp);
            
        else,%data length <=1 year
            
            sp=robustfit(time_array,data_array);
            
            v= sp(1) + sp(2) * time_array_min;
            
        end;
        
        
        inclination_secular_removed = inclination_est - v';
        
        
        % plot the data
        subplot(311);
        plot(time_array_min, f_data_secular_removed);
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        datetick('x','keeplimits');
        
        
        hold on;
        subplot(312);
        plot(time_array_min, declination_secular_removed);
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        datetick('x','keeplimits');
        
        hold on;
        subplot(313);
        plot(time_array_min, inclination_secular_removed);
        axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
        datetick('x','keeplimits');
        
        
        hold on;
        
    end;
    subplot(311);
    plot(time_array_min, scalar_f,'r');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    title(S(i).name(1:3));
    
    subplot(312);
    plot(time_array_min, declination_est_secular_removed,'r');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    subplot(313);
    plot(time_array_min, inclination_est_secular_removed,'r');
    axis([datenum(1995,1,1) datenum(2013,1,1) -inf inf]);
    datetick('x','keeplimits');
    
    
    
    %save the figure
    
    saveas(gcf,['/Users/manojnair/projects/obs_mag_data/Canada/plots/' S(i).name(1:3) '_F_DEC_DIP_compare'],'png');
    close all;
    clear mag_igrf DEC DIP H;
    
    % save the data
    
    save(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' ...
        S(i).name(1:3) '_F_DEC_DIP'], 'scalar_f', 'declination_est_secular_removed', ...
        'inclination_est_secular_removed', 'obj');
    
    clear mag_igrf DEC DIP H scalar_f declination declination_est declination_est_secular_removed ...
        declination_secular_removed x_data y_data z_data inclination inclination_est ...
        inclination_est_secular_removed inclination_secular_removed v time_array data_array b1;
    
    
    fprintf('Completed %s \n', S(i).name);
end;

%% Read the stations and calculate the RMS for X, Y and Z components



S = dir('/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/*.nc');

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);
load /Users/manojnair/projects/index/Ap_1995_2012.mat; % load daily Ap data
Ap_daily_rep_min = repmat(data(:,12),[1 1440]); % create repetition matrix
ap_data = reshape(Ap_daily_rep_min',[1,6575*1440]); % stretch the daily values to all the minutes in that dat
%LL = ap_data' < 0; % Kp > 6
LL = ap_data' < 20; % Exclude kp < 3


for i = 1:length(S),
    
    
    nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/' S(i).name] ;
    [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
    
    x_data_ref = x_data;
    y_data_ref = y_data;
    z_data_ref = z_data;
    
    
    L = isnan(x_data_ref) | LL ;
    
    if sum(~L) > 1000,
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , x_data_ref(~L)  );
        
        stats_x_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
        
    else,
        stats_x_ref(i,1:4) = NaN;
    end;
    
    L = isnan(y_data_ref) | LL;
    if sum(~L) > 1000,
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , y_data_ref(~L)  );
        
        stats_y_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
    else,
        stats_y_ref(i,1:4) = NaN;
    end;
    L = isnan(z_data_ref) | LL;
    if sum(~L) > 1000,
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , z_data_ref(~L)  );
        
        stats_z_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
    else,
        stats_z_ref(i,1:4) = NaN;
    end;
    
    
    for j = i+1:length(S),
        
        nc_fname = ['/Users/manojnair/projects/obs_mag_data/Canada/Secular_Removed_netCDF/' S(j).name] ;
        [x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);
        
        
        L = isnan(x_data) | isnan(x_data_ref)| LL;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , x_data(~L) - x_data_ref(~L) );
            
            stats_x(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_x(i,j,1:4) = NaN;
            
        end;
        
        L = isnan(y_data) | isnan(y_data_ref)| LL;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , y_data(~L) - y_data_ref(~L) );
            
            stats_y(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_y(i,j,1:4) = NaN;
            
        end;
        
        
        L = isnan(z_data) | isnan(z_data_ref)| LL;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , z_data(~L) - z_data_ref(~L) );
            
            stats_z(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_z(i,j,1:4) = NaN;
            
        end;
        
        fprintf('i=%2d j=%2d;', i, j);
    end;
    fprintf('Completed %s \n', S(i).name);
    
end;


%save /Users/manojnair/projects/correlation_length/NA_X_Y_Z_ApGT80_rms stats_x stats_x_ref stats_y stats_y_ref stats_z stats_z_ref;
save /Users/manojnair/projects/correlation_length/NA_X_Y_Z_ApGT20_rms stats_x stats_x_ref stats_y stats_y_ref stats_z stats_z_ref;

close all;
clear;

%% Read the stations and calculate the RMS of F between them



S =   dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*_F_DEC_DIP.mat')
time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
load /Users/manojnair/projects/index/Ap_1995_2012.mat; % load daily Ap data
Ap_daily_rep_min = repmat(data(:,12),[1 1440]); % create repetition matrix
ap_data = reshape(Ap_daily_rep_min',[1,6575*1440]); % stretch the daily values to all the minutes in that dat
%LL = ap_data' < 80; % Kp > 6
LL = ap_data' < 20; % Exclude kp < 3


for i = 1:length(S),
    
    
    load(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' S(i).name]);
    ref_scalar_f = scalar_f;
    ref_declination = declination_est_secular_removed;
    ref_inclination = inclination_est_secular_removed;
    
    
    L = isnan(ref_scalar_f) | LL ;
    
    if sum(~L) > 1000,
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , ref_scalar_f(~L)  );
        
        stats_f_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
    else
        stats_f_ref(i,1:4) = NaN;
    end;
    L = isnan(ref_declination)| LL ;
    if sum(~L) > 1000,
        
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , ref_declination(~L)  );
        
        stats_d_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
    else
        stats_d_ref(i,1:4) = NaN;
    end;
    L = isnan(ref_inclination)| LL ;
    if sum(~L) > 1000,
        
        
        [sp, stats_a] =    robustfit(time_array_min(~L) , ref_inclination(~L)  );
        
        stats_i_ref(i,:) = [stats_a.ols_s stats_a.robust_s stats_a.mad_s stats_a.s];
        
    else
        stats_i_ref(i,1:4) = NaN;
    end;
    
    for j = i+1:length(S),
        
        load(['/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/' S(j).name]);
        
        
        L = isnan(scalar_f) | isnan(ref_scalar_f)| LL ;;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , scalar_f(~L) - ref_scalar_f(~L) );
            
            stats_f(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_f(i,j,1:4) = NaN;
            
        end;
        
        L = isnan(declination_est_secular_removed) | isnan(ref_declination)| LL ;;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , declination_est_secular_removed(~L) - ref_declination(~L) );
            
            stats_d(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_d(i,j,1:4) = NaN;
            
        end;
        
        
        L = isnan(inclination_est_secular_removed) | isnan(ref_inclination)| LL ;;
        
        if sum(~L) > 1000,
            [sp, stats_diff] = robustfit(time_array_min(~L) , inclination_est_secular_removed(~L) - ref_inclination(~L) );
            
            stats_i(i,j,1:4) = [stats_diff.ols_s stats_diff.robust_s stats_diff.mad_s stats_diff.s];
            
        else,
            
            stats_i(i,j,1:4) = NaN;
            
        end;
        
        fprintf('i=%2d j=%2d;', i, j);
    end;
    fprintf('Completed %s \n', S(i).name);
    
end;


%save /Users/manojnair/projects/correlation_length/NA_F_DEC_DIP_ApGT80_rms stats_f stats_f_ref stats_d stats_d_ref stats_i stats_i_ref;
save /Users/manojnair/projects/correlation_length/NA_F_DEC_DIP_ApGT20_rms stats_f stats_f_ref stats_d stats_d_ref stats_i stats_i_ref;



%% find the distance in kilometers between stations

S =   dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*F_DEC_DIP.mat');
[master_a,master_b] = xlsread('/Users/manojnair/projects/obs_mag_data/Canada/Canada_array.xlsx');

master_b(1,:) = [];
station_list = cell2mat(master_b(:,1));


for i = 1:length(S),
    
    filestring(i,:) = S(i).name;
    
end;

[b, ia, ib] = intersect(station_list,S(i).name(1:3),'rows');

% Get the distance between stations

for i = 1:length(S) - 1,
    
    [b, ia, ib] = intersect(station_list,S(i).name(1:3),'rows');
    
    
    for j = i + 1 : length(S),
        
        [br, iar, ibr] = intersect(station_list,S(j).name(1:3),'rows');
        
        station_distance(i,j) = deg2km(distance(master_a(ia,1), ...
            master_a(ia,2),master_a(iar,1),master_a(iar,2)));
        
    end;
    fprintf('Completed %s \n', S(i).name);
    
end;

% Get Azimuths between stations

for i = 1:length(S),
    
    [b, ia, ib] = intersect(station_list,S(i).name(1:3),'rows');
    
    for j = 1 : length(S),
        
        [br, iar, ibr] = intersect(station_list,S(j).name(1:3),'rows');
        
        station_azimuth(i,j) =  azimuth(master_a(ia,1), ...
            master_a(ia,2),master_a(iar,1),master_a(iar,2),'degrees');
        
    end;
    fprintf('Completed %s \n', S(i).name);
    
end;




% Get the lat/long data of stations

for i = 1:length(S),
    
    [b, ia, ib] = intersect(station_list,S(i).name(1:3),'rows');
    
    station_lat(i)  = master_a(ia,1);
    station_lon(i) =  master_a(ia,2);
    station_gm_lat(i)  = master_a(ia,7);
    station_gm_lon(i) =  master_a(ia,8);
    station_iaga_code(i,:) = b;
    
    
end;

save /Users/manojnair/projects/correlation_length/NA_station_distance_SAH_removed  station_iaga_code station_azimuth station_distance station_lat station_lon station_gm_lat station_gm_lon;


%% Now plot data

S =   dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*_F_DEC_DIP.mat')

load /Users/manojnair/projects/correlation_length/NA_F_DEC_DIP_ApGT80_rms stats_i stats_i_ref;
load /Users/manojnair/projects/correlation_length/NA_X_Y_Z_ApGT20_rms stats_x stats_x_ref stats_y stats_y_ref stats_z stats_z_ref;
load /Users/manojnair/projects/correlation_length/NA_station_distance_SAH_removed station_distance;

stats = stats_i;
stats_ref = stats_i_ref;

nlength = length(S);

for i = 2:nlength-1,
    
    data = stats_z(1:i-1,i,2)';
    data(i:nlength-1) = stats_z(i,i+1:end,2);
    
    data1 = stats_x(1:i-1,i,2)';
    data1(i:nlength-1) = stats_x(i,i+1:end,2);
    
    
    sdist = station_distance(1:i-1,i)';
    sdist(i:nlength-1) = station_distance(i,i+1:end);
    normalized_rmse = data./sqrt(stats_z_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_z_ref(i,2).^2);
    normalized_rmse1 = data1./sqrt(stats_x_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_x_ref(i,2).^2);
    
    
    % discard data when station pairs are separated > 2000 km
    L = sdist > 2000;
    
    semilogx(sdist(~L),normalized_rmse(~L),'r.');
    hold on;
    semilogx(sdist(~L),normalized_rmse1(~L),'b.');
    hold off;
    
    
    %     L = isnan(normalized_rmse);
    %     [sp, stats_rms] =    robustfit(sdist(~L) , normalized_rmse(~L)  );
    %     hold on;
    %     loglog(sdist, sp(1) + sp(2) * sdist,'b.');
    %
    
    set(gca,'FontSize',16);
    %     text(50, 1, sprintf('%5.1f  %5.1f', sp(1), sp(2)));
    xlabel('log (Km)');
    ylabel('Normalized RMSE of DIP');
    axis([50  2100    0  2])
    title(S(i).name(1:3));
    %hold on;
    pause;
    
    close;
end;


%% plot maps

load /Users/manojnair/projects/correlation_length/NA_scalar_rms stats stats_ref;
load /Users/manojnair/projects/correlation_length/NA_station_distance station_distance  station_lat station_lon station_gm_lat station_gm_lon;
S =   dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*F.mat');

load coast;
nlength = length(S);

for i = 2:nlength-1,
    
    data = stats(1:i-1,i,2)';
    data(i:nlength-1) = stats(i,i+1:end,2);
    this_lat = station_lat([1:(i-1) i+1:nlength]);
    pair_lons = station_lon([1:(i-1) i+1:nlength]);
    normalized_rmse = data./sqrt(stats_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_ref(i,2).^2);
    sdist = station_distance(1:i-1,i)';
    sdist(i:nlength-1) = station_distance(i,i+1:end);
    
    subplot(122);
    
    worldmap([40,85],[-180,-50])
    plotm(lat,long);
    set(gca,'FontSize',16);
    
    
    for j = 1:length(normalized_rmse),
        
        if (~isnan(normalized_rmse(j)))
            
            plotm(this_lat(j),pair_lons(j),'r.', 'color', ...
                [[0 1 1]'*(1- normalized_rmse(j)./max(normalized_rmse))]','MarkerSize',30);
        end;
    end;
    plotm(station_lat(i),station_lon(i),'r.','MarkerSize',30);
    
    title(S(i).name(1:3));
    
    subplot(121);
    semilogx(sdist,normalized_rmse,'r.','MarkerSize',30);
    
    set(gca,'FontSize',16);
    xlabel('log (Km)');
    ylabel('Normalized RMSE of F');
    axis([10  10000    0.1  2])
    title(S(i).name(1:3));
    
    %hold on;
    pause;
    
    
end;



%% save the RMSE data (F, DIP , DEC, X, Y Z) to a file
% discard stations > 2000 km

load /Users/manojnair/projects/correlation_length/NA_station_distance_SAH_removed station_iaga_code station_distance station_azimuth station_lat station_lon station_gm_lat station_gm_lon;
S =   dir('/Users/manojnair/projects/obs_mag_data/Canada/work/station_matlab_files_baseline_removed/*F_DEC_DIP.mat');
load /Users/manojnair/projects/correlation_length/NA_F_DEC_DIP_ApGT20_rms.mat stats_f stats_f_ref stats_d stats_d_ref stats_i stats_i_ref;
load /Users/manojnair/projects/correlation_length/NA_X_Y_Z_ApGT20_rms stats_x stats_x_ref stats_y stats_y_ref stats_z stats_z_ref;

fid = fopen('/Users/manojnair/projects/correlation_length/North_America_rmse.txt','wt');

fid_station_list = fopen('/Users/manojnair/projects/correlation_length/North_America_list.txt','wt');

fprintf(fid,['ID1, LAT1, LON1, ID2, LAT2, LON2, DISTANCE(Km), AZIMUTH(degrees), RMSE_F, RMS_F_1, RMS_F_2,' ...
    ' RMSE_DEC, RMS_DEC_1, RMS_DEC_2, RMSE_DIP, RMS_DIP_1, RMS_DIP_2, RMSE_X, RMS_X_1, RMS_X_2,' ...
    ' RMSE_Y, RMS_Y_1, RMS_Y_2, RMSE_Z RMS_Z_1, RMS_Z_2,\n']);

nlength = length(S);

for i = 2:nlength-1,
%for i = [2:41 43:62],
    
    % Get the lat/long of base station
    this_lat = station_lat(i);
    this_lon = station_lon(i);
    
    % get rms of different components for this station
    this_f = stats_f_ref(i,2);
    this_x = stats_x_ref(i,2);
    this_y = stats_y_ref(i,2);
    this_z = stats_z_ref(i,2);
    this_d = stats_d_ref(i,2);
    this_i = stats_i_ref(i,2);
    
    % get rms of the paired stations
    paired_f = stats_f_ref([1:(i-1) i+1:nlength],2);
    paired_x = stats_x_ref([1:(i-1) i+1:nlength],2);
    paired_y = stats_y_ref([1:(i-1) i+1:nlength],2);
    paired_z = stats_z_ref([1:(i-1) i+1:nlength],2);
    paired_d = stats_d_ref([1:(i-1) i+1:nlength],2);
    paired_i = stats_i_ref([1:(i-1) i+1:nlength],2);
    
    
    
    % Get the location of paired stations
    paired_lats = station_lat([1:(i-1) i+1:nlength]);
    paired_lons = station_lon([1:(i-1) i+1:nlength]);
    
    %Get the IAGA code for this station
    this_code = station_iaga_code(i,:);
    
    %Get IAGA codes for paired stations
    paired_codes = station_iaga_code([1:(i-1) i+1:nlength],:);
    % Get the distance to all stations
    sdist = station_distance(1:i-1,i)';
    sdist(i:nlength-1) = station_distance(i,i+1:end);
    % Get the azimuth of remote station for each pair
    sazimuth = station_azimuth(i,[1:(i-1) i+1:nlength]);
    
    
    % Calculate normalized RMSE of F
    data = stats_f(1:i-1,i,2)';
    data(i:nlength-1) = stats_f(i,i+1:end,2)';
    %normalized_f_rmse = data./sqrt(stats_f_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_f_ref(i,2).^2);
    f_rmse = data;
    
    % Calculate normalized rmse of decliantiom
    
    data = stats_d(1:i-1,i,2)';
    data(i:nlength-1) = stats_d(i,i+1:end,2);
    %normalized_d_rmse = data./sqrt(stats_d_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_d_ref(i,2).^2);
    d_rmse = data;
    
    % calculate normalized rmse of inclination
    
    data = stats_i(1:i-1,i,2)';
    data(i:nlength-1) = stats_i(i,i+1:end,2);
    %normalized_i_rmse = data./sqrt(stats_i_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_i_ref(i,2).^2);
    i_rmse = data;
    
    % calculate normalized rmse of X
    
    data = stats_x(1:i-1,i,2)';
    data(i:nlength-1) = stats_x(i,i+1:end,2);
    %normalized_x_rmse = data./sqrt(stats_x_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_x_ref(i,2).^2);
    x_rmse = data;
    
    % calculate normalized rmse of Y
    
    data = stats_y(1:i-1,i,2)';
    data(i:nlength-1) = stats_y(i,i+1:end,2);
    %normalized_y_rmse = data./sqrt(stats_y_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_y_ref(i,2).^2);
    y_rmse = data;
    
    % calculate normalized rmse of Z
    
    data = stats_z(1:i-1,i,2)';
    data(i:nlength-1) = stats_z(i,i+1:end,2);
    %normalized_z_rmse = data./sqrt(stats_z_ref([1:(i-1) i+1:nlength],2)'.^2 + stats_z_ref(i,2).^2);
    z_rmse = data;
    
    
    % Discard data for station pairs located at > 2000 km | no data
    
    L = sdist > 2000 | isnan(f_rmse) | isnan(d_rmse) ...
        | isnan(i_rmse) | isnan(z_rmse) | isnan(y_rmse)  ...
        | isnan(x_rmse);
    
    paired_codes(L,:) = [];
    paired_lats(L) = [];
    paired_lons(L) = [];
    sdist(L) = [];
    sazimuth(L) = [];
    %     normalized_f_rmse(L) = [];
    %     normalized_d_rmse(L) = [];
    %     normalized_i_rmse(L) = [];
    %     normalized_x_rmse(L) = [];
    %     normalized_y_rmse(L) = [];
    %     normalized_z_rmse(L) = [];
    f_rmse(L) = [];
    x_rmse(L) = [];
    y_rmse(L) = [];
    z_rmse(L) = [];
    d_rmse(L) = [];
    i_rmse(L) = [];
    paired_f(L) = [];
    paired_x(L) = [];
    paired_y(L) = [];
    paired_z(L) = [];
    paired_d(L) = [];
    paired_i(L) = [];
    
    
    % Now, print the values ID1, LAT1,LON1, ID2, LAT2, LON2, DISTANCE, AZIMUTH,
    % F, DEC, DIP, X, Y, Z
    
    %   fprintf(fid,['ID1, LAT1,LON1, ID2, LAT2, LON2, DISTANCE(Km), AZIMUTH(degrees), NRMSE_F, RMS_F_1, RMS_F_2,' ...
    % 'NRMSE_DEC, RMS_DEC_1, RMS_DEC_2, NRMSE_DIP, RMS_DIP_1, RMS_DIP_2, NRMSE_X, RMS_X_1, RMS_X_2,' ...
    % 'NRMSE_Y, RMS_Y_1, RMS_Y_2, NRMSE_Z RMS_Z_1, RMS_Z_2,\n']);
    
    
    
    for j = 1:length(sdist),
        
        fprintf(fid,'%s %5.2f %5.2f %s %5.2f %5.2f %6.1f %5.2f %5.3f %5.3f %5.3f %10.7f %10.7f %10.7f %10.7f %10.7f %10.7f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f %5.3f\n' ...
            , this_code, this_lat, this_lon, paired_codes(j,:), paired_lats(j), ...
            paired_lons(j), sdist(j), sazimuth(j), f_rmse(j), this_f, paired_f(j), ...
            d_rmse(j), this_d, paired_d(j),  i_rmse(j), this_i, paired_i(j), ...
            x_rmse(j), this_x, paired_x(j),  y_rmse(j), this_y, paired_y(j), ...
            z_rmse(j), this_z, paired_z(j));
        
    end;
    
    fprintf(fid_station_list,'%s %5.2f %5.2f\n', this_code, this_lat, this_lon);
    
    clear data this_code this_lat this_lon paired_codes paired_lats ...
        paired_lons sdist sazimuth normalized_f_rmse ...
        normalized_d_rmse normalized_i_rmse normalized_x_rmse ...
        normalized_y_rmse normalized_z_rmse;
end;

fclose(fid);
fclose(fid_station_list);



